JS23.2_QC
Jaime So
2023-10-04
Perform Library Quality Control Checks
Before running enrichment analysis we need to make sure each sample meets read depth and alignment quality cutoffs.
I wrote three functions to perform preliminary plate QC.
- bt_align() merges bowtie standard error report files with sample well information for each plate
- plot_96well() takes bt_align output df and plots an overview of total read depth and % alignment for each sample on a plate
- import.counts() reads in htseq counts files and makes a read count
matrix while applying quality depth and alignment cutoffs.
- samples with <10*(#represented phage) reads or <70% reads aligned are automatically thrown out
- additionally, mock samples that do not meet specified pearson’s correlation cutoff (<0.8) are excluded from the matrix
Generate final counts matrices for each plate for input into edgeR enrichment analysis
Read botwie alignment standard error reports
p1 <- bt_align("p1_PhIP_out/bt_out/", "../Plate_spreadsheets/plate1_wells.xlsx")
p2 <- bt_align("p2-3_PhIP_out/bt_out/", "../Plate_spreadsheets/plate2_wells.xlsx")
p3 <- bt_align("p2-3_PhIP_out/bt_out/", "../Plate_spreadsheets/plate3_wells.xlsx")
p4 <- bt_align("p4-5_PhIP_out/bt_out/", "../Plate_spreadsheets/plate4_wells.xlsx")
p5 <- bt_align("p4-5_PhIP_out/bt_out/", "../Plate_spreadsheets/plate5_wells.xlsx")
p6 <- bt_align("PhIP_out_Tbg6-8_Tbr9/bt_out/", "../Plate_spreadsheets/plate6_wells.xlsx")
p7 <- bt_align("PhIP_out_Tbg6-8_Tbr9/bt_out/", "../Plate_spreadsheets/plate7_wells.xlsx")
p8 <- bt_align("PhIP_out_Tbg6-8_Tbr9/bt_out/", "../Plate_spreadsheets/plate8_wells.xlsx")Plate 1
plot_96well(p1)
plate1 <- import.counts(directory = "p1_PhIP_out/htseq_out/",
"_L002_R1_001_trimmed.fq.btSE.txt", p1, 0.8)
## b_1 b_2 b_3 b_4 b_5 b_6 b_7 b_8
## beads_r1 1
## beads_r2 , 1
## beads_r3 + , 1
## beads_r4 * , * 1
## beads_r5 * , + B 1
## beads_r6 * , * B * 1
## beads_r7 + , + + + * 1
## beads_r8 * , + B * * * 1
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1Plate 2
plot_96well(p2)
plate2 <- import.counts(directory = "p2-3_PhIP_out/",
"_L002_R1_001_trimmed.fq.btSE.txt", p2, 0.8)
## b_2_1 b_2_2 b_2_3 b_2_4 b_2_5 b_2_6 b_2_7 b_2_8
## beads_p2_r1 1
## beads_p2_r2 * 1
## beads_p2_r3 B * 1
## beads_p2_r4 B * B 1
## beads_p2_r5 B * B B 1
## beads_p2_r6 B * B B B 1
## beads_p2_r7 B * B B B B 1
## beads_p2_r8 B * B B B B B 1
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1Plate 3
plot_96well(p3)
plate3 <- import.counts(directory = "p2-3_PhIP_out/",
"_L002_R1_001_trimmed.fq.btSE.txt", p3, 0.8)
## b_3_1 b_3_2 b_3_3 b_3_4 b_3_5 b_3_6 b_3_7 b_3_8
## beads_p3_r1 1
## beads_p3_r2 + 1
## beads_p3_r3 B + 1
## beads_p3_r4 * + B 1
## beads_p3_r5 + , + + 1
## beads_p3_r6 B + B B + 1
## beads_p3_r7 B + B * + B 1
## beads_p3_r8 , , , , , , , 1
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1Plate 4
plot_96well(p4)
plate4 <- import.counts(directory = "p4-5_PhIP_out/htseq_out/",
"_L002_R1_001_trimmed.fq.btSE.txt", p4, 0.8)
## b_4_1 b_4_2 b_4_3 b_4_4 b_4_5 b_4_6 b_4_7 b_4_8
## beads_p4_r1 1
## beads_p4_r2 * 1
## beads_p4_r3 * * 1
## beads_p4_r4 * * * 1
## beads_p4_r5 * * B * 1
## beads_p4_r6 * + + + + 1
## beads_p4_r7 * + * * * * 1
## beads_p4_r8 * * B * * + * 1
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1Plate 5
plot_96well(p5)
plate5 <- import.counts(directory = "p4-5_PhIP_out/htseq_out/",
"_L002_R1_001_trimmed.fq.btSE.txt", p5, 0.8)
## b_5_1 b_5_2 b_5_3 b_5_4 b_5_5 b_5_6 b_5_7 b_5_8
## beads_p5_r1 1
## beads_p5_r2 , 1
## beads_p5_r3 , , 1
## beads_p5_r4 , , + 1
## beads_p5_r5 + + + + 1
## beads_p5_r6 + , + + * 1
## beads_p5_r7 + , + + * * 1
## beads_p5_r8 + , + + * * * 1
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1Plate 6
plot_96well(p6)
plate6 <- import.counts(directory = "PhIP_out_Tbg6-8_Tbr9/htseq_out/lane1/",
"_L001_R1_001_trimmed.fq.btSE.txt", p6, 0.8)
## b_6_1 b_6_2 b_6_3 b_6_4 b_6_5 b_6_6 b_6_7 b_6_8
## beads_p6_r1 1
## beads_p6_r2 + 1
## beads_p6_r3 + + 1
## beads_p6_r4 + + * 1
## beads_p6_r5 + + * * 1
## beads_p6_r6 , + , , , 1
## beads_p6_r7 + + * * * , 1
## beads_p6_r8 + + * * * , * 1
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1Plate 7
plot_96well(p7)
plate7 <- import.counts(directory = "PhIP_out_Tbg6-8_Tbr9/htseq_out/lane1/",
"_L001_R1_001_trimmed.fq.btSE.txt", p7, 0.8)
## b_7_1 b_7_2 b_7_3 b_7_4 b_7_5 b_7_6 b_7_7 b_7_8
## beads_p7_r1 1
## beads_p7_r2 * 1
## beads_p7_r3 * B 1
## beads_p7_r4 , , , 1
## beads_p7_r5 + + + + 1
## beads_p7_r6 * B B , + 1
## beads_p7_r7 * B B , + B 1
## beads_p7_r8 , , , * , , , 1
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1Plate 8
plot_96well(p8)
plate8 <- import.counts(directory = "PhIP_out_Tbg6-8_Tbr9/htseq_out/lane2/",
"_L002_R1_001_trimmed.fq.btSE.txt", p8, 0.8)
## b_8_1 b_8_2 b_8_3 b_8_4 b_8_5 b_8_6 b_8_7 b_8_8
## beads_p8_r1 1
## beads_p8_r2 B 1
## beads_p8_r3 B B 1
## beads_p8_r4 B B B 1
## beads_p8_r5 B B B B 1
## beads_p8_r6 B B B B B 1
## beads_p8_r7 B B B B B B 1
## beads_p8_r8 , , , , , , , 1
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1Save the Rdata space